# Replication Archive for: 
# Coppock, Alexander and Donald P. Green. 2020. 
# "Do Belief Systems Exhibit Dynamic Constraint?" 
# The Journal of Politics, Forthcoming.

rm(list = ls())

library(tidyverse)
library(car)
library(psych)

# Load deduped, anonymized datasets

study_1_mturk_w1 <- read_csv("data/raw/study_1_mturk_wave_1_raw.csv")
study_1_mturk_w2 <- read_csv("data/raw/study_1_mturk_wave_2_raw.csv")
study_1_mturk_w3 <- read_csv("data/raw/study_1_mturk_wave_3_raw.csv")

# Wave 1 cleaning ---------------------------------------------------------

study_1_mturk_w1 <- within(study_1_mturk_w1, {
  
  # Covariates
  
  rog_type <- rep(NA, nrow(study_1_mturk_w1))
  rog_type[ROG_1 == 1 & ROG_2 == 1] <- "Conservative"
  rog_type[ROG_1 == 1 & ROG_2 == 2] <- "Libertarian"
  rog_type[ROG_1 == 2 & ROG_2 == 1] <- "Communitarian"
  rog_type[ROG_1 == 2 & ROG_2 == 2] <- "Liberal"
  
  pid_7 <- rep(NA, nrow(study_1_mturk_w1))
  pid_7[pid_3 == 1 & pid_dem == 1] <- 1
  pid_7[pid_3 == 1 & pid_dem == 2] <- 2
  pid_7[pid_3 == 1 & is.na(pid_dem)] <- 2
  pid_7[pid_3 == 1 & pid_dem == ""] <- 2
  
  pid_7[pid_3 == 4] <- 4
  
  pid_7[pid_3 == 2 & pid_ind == 3] <- 3
  pid_7[pid_3 == 2 & pid_ind == 4] <- 4
  pid_7[pid_3 == 2 & is.na(pid_ind)] <- 4
  pid_7[pid_3 == 2 & pid_ind == 5] <- 5
  
  pid_7[pid_3 == 3 & is.na(pid_rep)] <- 6
  pid_7[pid_3 == 3 & pid_rep == 6] <- 6
  pid_7[pid_3 == 3 & pid_rep == 7] <- 7
  
  pid_3_cat <- recode(pid_7, "1:3='Democrat';4='Independent';5:7='Republican'")
  
  ideo_5 <- recode(ideology, "1='Liberal';2='Moderate';3='Conservative';4='Libertarian';5:6='Other'")
  ideo_5 <- factor(ideo_5, levels = c("Liberal", "Moderate", "Libertarian", "Conservative", "Other"))
  
  educ_5 <- as.numeric(education)
  educ_5[is.na(education)] <- 3
  
  educ_5 <- factor(educ_5, levels = 1:5, labels = c("Less than High School", "High School", "Some College", "College", "Graduate School"))
  
  race_4 <- rep(NA, nrow(study_1_mturk_w1))
  race_4[race == 1] <- "White"
  race_4[race == 2] <- "Black"
  race_4[race == 3] <- "Hispanic"
  race_4[race %in% 4:7] <- "Other"
  
  race_4 <- factor(race_4, levels = c("Black", "Hispanic", "White", "Other"))
  
  female <- ifelse(gender==2, 1, 0)
  
  age_5 <- as.numeric(age)
  age_5[age==9] <-4
  age_5[age==8] <-7
  age_5 <- age_5 - 2
  
  age_5 <- factor(age_5, levels = 1:5, labels = c("18 - 29", "30 - 39", "40 - 49", "50 - 59", "60+"))
  
  
  Z <- relevel(factor(Z), ref="control")
  
  #Recode Outcome Measures so larger # indicate more agreement with author
  
  dv_climate_1_w1 <- as.numeric(climate_1_1)
  dv_climate_2_w1 <- as.numeric(climate_2_1)
  dv_climate_3_w1 <- -1*as.numeric(climate_3_1) + 8
  dv_climate_4_w1 <- -1*as.numeric(climate_4_1) + 8
  dv_climate_5_w1 <- -1*as.numeric(climate_5_1) 
  
  dv_wall_1_w1 <- as.numeric(wall_1_1)
  dv_wall_2_w1 <- -1*as.numeric(wall_2_1) + 101
  dv_wall_3_w1 <- -1*as.numeric(wall_3_1) + 8
  dv_wall_4_w1 <- as.numeric(wall_4_1)
  
  dv_amtrak_1_w1 <- as.numeric(amtrak_1_1)
  dv_amtrak_1_w1[as.numeric(amtrak_1_1) %in% c(6,7,8)] <- 
    dv_amtrak_1_w1[as.numeric(amtrak_1_1) %in% c(6,7,8)] -1
  # hard brackets subsetting, fixing it since the raw data is missing a 5
  dv_amtrak_2_w1 <- as.numeric(amtrak_2_1)
  dv_amtrak_3_w1 <- -1*as.numeric(amtrak_3_1) + 8
  dv_amtrak_4_w1 <- as.numeric(amtrak_4_1)
  
  dv_vets_1_w1 <- -1*as.numeric(vets_1_1)+101
  dv_vets_2_w1 <- as.numeric(vets_2_1)
  dv_vets_3_w1 <- -1*as.numeric(vets_3_1) + 8
  dv_vets_4_w1 <- as.numeric(vets_4_1)
  
  dv_flat_1_w1 <- -1 * as.numeric(flat_1_1) + 8
  dv_flat_2_w1 <- -1 * as.numeric(flat_2_1) + 101
  dv_flat_3_w1 <- -1 * as.numeric(flat_3_1) + 8
  dv_flat_4_w1 <- -1 * as.numeric(flat_4_1) + 8
  
  dv_climate_main_w1 <- dv_climate_1_w1
  dv_wall_main_w1 <- dv_wall_3_w1
  dv_amtrak_main_w1 <- dv_amtrak_1_w1
  dv_vets_main_w1 <- dv_vets_2_w1
  dv_flat_main_w1 <- dv_flat_1_w1
  
})

study_1_mturk_w1 <-
  study_1_mturk_w1 %>%
  select(
    subject_id,
    time_spent_reading,
    rog_type,
    pid_7,
    pid_3_cat,
    ideo_5,
    educ_5,
    race_4,
    female,
    age_5,
    Z,
    starts_with("dv_")
  )

# Wave 2 cleaning ---------------------------------------------------------

study_1_mturk_w2 <- within(study_1_mturk_w2,{
  
  dv_climate_1_w2 <- as.numeric(climate_1_1)
  dv_climate_2_w2 <- as.numeric(climate_2_1)
  dv_climate_3_w2 <- -1*as.numeric(climate_3_1) + 8
  dv_climate_4_w2 <- -1*as.numeric(climate_4_1) + 8
  dv_climate_5_w2 <- -1*as.numeric(climate_5_1) 
  
  dv_wall_1_w2 <- as.numeric(wall_1_1)
  dv_wall_2_w2 <- -1*as.numeric(wall_2_1) + 101
  dv_wall_3_w2 <- -1*as.numeric(wall_3_1) + 8
  dv_wall_4_w2 <- as.numeric(wall_4_1)
  
  dv_amtrak_1_w2 <- as.numeric(amtrak_1_1)
  dv_amtrak_1_w2[as.numeric(amtrak_1_1) %in% c(6,7,8)] <- 
    dv_amtrak_1_w2[as.numeric(amtrak_1_1) %in% c(6,7,8)] -1
  dv_amtrak_2_w2 <- as.numeric(amtrak_2_1)
  dv_amtrak_3_w2 <- -1*as.numeric(amtrak_3_1) + 8
  dv_amtrak_4_w2 <- as.numeric(amtrak_4_1)
  
  dv_vets_1_w2 <- -1*as.numeric(vets_1_1)+101
  dv_vets_2_w2 <- as.numeric(vets_2_1)
  dv_vets_3_w2 <- -1*as.numeric(vets_3_1) + 8
  dv_vets_4_w2 <- as.numeric(vets_4_1)
  
  dv_flat_1_w2 <- -1 * as.numeric(flat_1_1) + 8
  dv_flat_2_w2 <- -1 * as.numeric(flat_2_1) +101
  dv_flat_3_w2 <- -1* as.numeric(flat_3_1) + 8
  dv_flat_4_w2 <- -1 * as.numeric(flat_4_1) + 8
  
  dv_climate_main_w2 <- dv_climate_1_w2
  dv_wall_main_w2 <- dv_wall_3_w2
  dv_amtrak_main_w2 <- dv_amtrak_1_w2
  dv_vets_main_w2 <- dv_vets_2_w2
  dv_flat_main_w2 <- dv_flat_1_w2
  
  # Distraction Randomization
  Z_distract <- as.numeric(Z_distract)
  
  dv_mj_1_w2 <- -1 * as.numeric(mj_1_1) + 8
  dv_mj_2_w2 <- -1 * as.numeric(mj_2_1) + 8
  dv_mj_3_w2 <- -1 * as.numeric(mj_3_2) + 8
  dv_mj_4_w2 <- -1 * as.numeric(mj_4_2) + 8
  
})

study_1_mturk_w2 <-
  study_1_mturk_w2 %>%
  select(subject_id,
         Z_distract,
         starts_with("dv_"))

# Wave 3 cleaning ---------------------------------------------------------

study_1_mturk_w3 <- within(study_1_mturk_w3, {
  dv_climate_1_w3 <- as.numeric(climate_1_1)
  dv_climate_2_w3 <- as.numeric(climate_2_1)
  dv_climate_3_w3 <- -1 * as.numeric(climate_3_1) + 8
  dv_climate_4_w3 <- -1 * as.numeric(climate_4_1) + 8
  dv_climate_5_w3 <- -1 * as.numeric(climate_5_1)
  
  dv_wall_1_w3 <- as.numeric(wall_1_1)
  dv_wall_2_w3 <- -1 * as.numeric(wall_2_1) + 101
  dv_wall_3_w3 <- -1 * as.numeric(wall_3_1) + 8
  dv_wall_4_w3 <- as.numeric(wall_4_1)
  
  dv_amtrak_1_w3 <- as.numeric(amtrak_1_1)
  dv_amtrak_1_w3[as.numeric(amtrak_1_1) %in% c(6, 7, 8)] <-
    dv_amtrak_1_w3[as.numeric(amtrak_1_1) %in% c(6, 7, 8)] - 1
  dv_amtrak_2_w3 <- as.numeric(amtrak_2_1)
  dv_amtrak_3_w3 <- -1 * as.numeric(amtrak_3_1) + 8
  dv_amtrak_4_w3 <- as.numeric(amtrak_4_1)
  
  dv_vets_1_w3 <- -1 * as.numeric(vets_1_1) + 101
  dv_vets_2_w3 <- as.numeric(vets_2_1)
  dv_vets_3_w3 <- -1 * as.numeric(vets_3_1) + 8
  dv_vets_4_w3 <- as.numeric(vets_4_1)
  
  dv_flat_1_w3 <- -1 * as.numeric(flat_1_1) + 8
  dv_flat_2_w3 <- -1 * as.numeric(flat_2_1) + 101
  dv_flat_3_w3 <- -1 * as.numeric(flat_3_1) + 8
  dv_flat_4_w3 <- -1 * as.numeric(flat_4_1) + 8
  
  dv_climate_main_w3 <- dv_climate_1_w3
  dv_wall_main_w3 <- dv_wall_3_w3
  dv_amtrak_main_w3 <- dv_amtrak_1_w3
  dv_vets_main_w3 <- dv_vets_2_w3
  dv_flat_main_w3 <- dv_flat_1_w3
  
  dv_mj_1_w3 <- -1 * as.numeric(mj_1_1) + 8
  dv_mj_2_w3 <- -1 * as.numeric(mj_2_1) + 8
  dv_mj_3_w3 <- -1 * as.numeric(mj_3_2) + 8
  dv_mj_4_w3 <- -1 * as.numeric(mj_4_2) + 8
})

study_1_mturk_w3 <-
  study_1_mturk_w3 %>%
  select(subject_id,
         starts_with("dv_"))


# Merge datasets ----------------------------------------------------------

study_1_mturk <-
  study_1_mturk_w1 %>%
  left_join(study_1_mturk_w2) %>%
  left_join(study_1_mturk_w3)



# Standardize Outcomes ----------------------------------------------------


study_1_mturk <- within(study_1_mturk,{
  
  dv_climate_1_s_w1 <- dv_climate_1_w1/sd(dv_climate_1_w1[Z=="control"],na.rm = TRUE)
  dv_climate_2_s_w1 <- dv_climate_2_w1/sd(dv_climate_2_w1[Z=="control"],na.rm = TRUE)
  dv_climate_3_s_w1 <- dv_climate_3_w1/sd(dv_climate_3_w1[Z=="control"],na.rm = TRUE)
  dv_climate_4_s_w1 <- dv_climate_4_w1/sd(dv_climate_4_w1[Z=="control"],na.rm = TRUE)
  dv_climate_5_s_w1 <- dv_climate_5_w1/sd(dv_climate_5_w1[Z=="control"],na.rm = TRUE)
  dv_wall_1_s_w1 <- dv_wall_1_w1/sd(dv_wall_1_w1[Z=="control"],na.rm = TRUE)
  dv_wall_2_s_w1 <- dv_wall_2_w1/sd(dv_wall_2_w1[Z=="control"],na.rm = TRUE)
  dv_wall_3_s_w1 <- dv_wall_3_w1/sd(dv_wall_3_w1[Z=="control"],na.rm = TRUE)
  dv_wall_4_s_w1 <- dv_wall_4_w1/sd(dv_wall_4_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_1_s_w1 <- dv_amtrak_1_w1/sd(dv_amtrak_1_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_2_s_w1 <- dv_amtrak_2_w1/sd(dv_amtrak_2_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_3_s_w1 <- dv_amtrak_3_w1/sd(dv_amtrak_3_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_4_s_w1 <- dv_amtrak_4_w1/sd(dv_amtrak_4_w1[Z=="control"],na.rm = TRUE)
  dv_vets_1_s_w1 <- dv_vets_1_w1/sd(dv_vets_1_w1[Z=="control"],na.rm = TRUE)
  dv_vets_2_s_w1 <- dv_vets_2_w1/sd(dv_vets_2_w1[Z=="control"],na.rm = TRUE)
  dv_vets_3_s_w1 <- dv_vets_3_w1/sd(dv_vets_3_w1[Z=="control"],na.rm = TRUE)
  dv_vets_4_s_w1 <- dv_vets_4_w1/sd(dv_vets_4_w1[Z=="control"],na.rm = TRUE)
  dv_flat_1_s_w1 <- dv_flat_1_w1/sd(dv_flat_1_w1[Z=="control"],na.rm = TRUE)
  dv_flat_2_s_w1 <- dv_flat_2_w1/sd(dv_flat_2_w1[Z=="control"],na.rm = TRUE)
  dv_flat_3_s_w1 <- dv_flat_3_w1/sd(dv_flat_3_w1[Z=="control"],na.rm = TRUE)
  dv_flat_4_s_w1 <- dv_flat_4_w1/sd(dv_flat_4_w1[Z=="control"],na.rm = TRUE)
  dv_climate_main_s_w1 <- dv_climate_1_s_w1
  dv_wall_main_s_w1 <- dv_wall_3_s_w1
  dv_amtrak_main_s_w1 <- dv_amtrak_1_s_w1
  dv_vets_main_s_w1 <- dv_vets_2_s_w1
  dv_flat_main_s_w1 <- dv_flat_1_s_w1
  
  dv_climate_1_s_w2 <- dv_climate_1_w2/sd(dv_climate_1_w1[Z=="control"],na.rm = TRUE)
  dv_climate_2_s_w2 <- dv_climate_2_w2/sd(dv_climate_2_w1[Z=="control"],na.rm = TRUE)
  dv_climate_3_s_w2 <- dv_climate_3_w2/sd(dv_climate_3_w1[Z=="control"],na.rm = TRUE)
  dv_climate_4_s_w2 <- dv_climate_4_w2/sd(dv_climate_4_w1[Z=="control"],na.rm = TRUE)
  dv_climate_5_s_w2 <- dv_climate_5_w2/sd(dv_climate_5_w1[Z=="control"],na.rm = TRUE)
  dv_wall_1_s_w2 <- dv_wall_1_w2/sd(dv_wall_1_w1[Z=="control"],na.rm = TRUE)
  dv_wall_2_s_w2 <- dv_wall_2_w2/sd(dv_wall_2_w1[Z=="control"],na.rm = TRUE)
  dv_wall_3_s_w2 <- dv_wall_3_w2/sd(dv_wall_3_w1[Z=="control"],na.rm = TRUE)
  dv_wall_4_s_w2 <- dv_wall_4_w2/sd(dv_wall_4_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_1_s_w2 <- dv_amtrak_1_w2/sd(dv_amtrak_1_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_2_s_w2 <- dv_amtrak_2_w2/sd(dv_amtrak_2_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_3_s_w2 <- dv_amtrak_3_w2/sd(dv_amtrak_3_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_4_s_w2 <- dv_amtrak_4_w2/sd(dv_amtrak_4_w1[Z=="control"],na.rm = TRUE)
  dv_vets_1_s_w2 <- dv_vets_1_w2/sd(dv_vets_1_w1[Z=="control"],na.rm = TRUE)
  dv_vets_2_s_w2 <- dv_vets_2_w2/sd(dv_vets_2_w1[Z=="control"],na.rm = TRUE)
  dv_vets_3_s_w2 <- dv_vets_3_w2/sd(dv_vets_3_w1[Z=="control"],na.rm = TRUE)
  dv_vets_4_s_w2 <- dv_vets_4_w2/sd(dv_vets_4_w1[Z=="control"],na.rm = TRUE)
  dv_flat_1_s_w2 <- dv_flat_1_w2/sd(dv_flat_1_w1[Z=="control"],na.rm = TRUE)
  dv_flat_2_s_w2 <- dv_flat_2_w2/sd(dv_flat_2_w1[Z=="control"],na.rm = TRUE)
  dv_flat_3_s_w2 <- dv_flat_3_w2/sd(dv_flat_3_w1[Z=="control"],na.rm = TRUE)
  dv_flat_4_s_w2 <- dv_flat_4_w2/sd(dv_flat_4_w1[Z=="control"],na.rm = TRUE)

  dv_mj_1_s_w2 <- dv_mj_1_w2/sd(dv_mj_1_w2[Z_distract==0],na.rm = TRUE)
  dv_mj_2_s_w2 <- dv_mj_2_w2/sd(dv_mj_2_w2[Z_distract==0],na.rm = TRUE)
  dv_mj_3_s_w2 <- dv_mj_3_w2/sd(dv_mj_3_w2[Z_distract==0],na.rm = TRUE)
  dv_mj_4_s_w2 <- dv_mj_4_w2/sd(dv_mj_4_w2[Z_distract==0],na.rm = TRUE)
  
  dv_climate_main_s_w2 <- dv_climate_1_s_w2
  dv_wall_main_s_w2 <- dv_wall_3_s_w2
  dv_amtrak_main_s_w2 <- dv_amtrak_1_s_w2
  dv_vets_main_s_w2 <- dv_vets_2_s_w2
  dv_flat_main_s_w2 <- dv_flat_1_s_w2
  
  dv_climate_1_s_w3 <- dv_climate_1_w3/sd(dv_climate_1_w1[Z=="control"],na.rm = TRUE)
  dv_climate_2_s_w3 <- dv_climate_2_w3/sd(dv_climate_2_w1[Z=="control"],na.rm = TRUE)
  dv_climate_3_s_w3 <- dv_climate_3_w3/sd(dv_climate_3_w1[Z=="control"],na.rm = TRUE)
  dv_climate_4_s_w3 <- dv_climate_4_w3/sd(dv_climate_4_w1[Z=="control"],na.rm = TRUE)
  dv_climate_5_s_w3 <- dv_climate_5_w3/sd(dv_climate_5_w1[Z=="control"],na.rm = TRUE)
  dv_wall_1_s_w3 <- dv_wall_1_w3/sd(dv_wall_1_w1[Z=="control"],na.rm = TRUE)
  dv_wall_2_s_w3 <- dv_wall_2_w3/sd(dv_wall_2_w1[Z=="control"],na.rm = TRUE)
  dv_wall_3_s_w3 <- dv_wall_3_w3/sd(dv_wall_3_w1[Z=="control"],na.rm = TRUE)
  dv_wall_4_s_w3 <- dv_wall_4_w3/sd(dv_wall_4_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_1_s_w3 <- dv_amtrak_1_w3/sd(dv_amtrak_1_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_2_s_w3 <- dv_amtrak_2_w3/sd(dv_amtrak_2_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_3_s_w3 <- dv_amtrak_3_w3/sd(dv_amtrak_3_w1[Z=="control"],na.rm = TRUE)
  dv_amtrak_4_s_w3 <- dv_amtrak_4_w3/sd(dv_amtrak_4_w1[Z=="control"],na.rm = TRUE)
  dv_vets_1_s_w3 <- dv_vets_1_w3/sd(dv_vets_1_w1[Z=="control"],na.rm = TRUE)
  dv_vets_2_s_w3 <- dv_vets_2_w3/sd(dv_vets_2_w1[Z=="control"],na.rm = TRUE)
  dv_vets_3_s_w3 <- dv_vets_3_w3/sd(dv_vets_3_w1[Z=="control"],na.rm = TRUE)
  dv_vets_4_s_w3 <- dv_vets_4_w3/sd(dv_vets_4_w1[Z=="control"],na.rm = TRUE)
  dv_flat_1_s_w3 <- dv_flat_1_w3/sd(dv_flat_1_w1[Z=="control"],na.rm = TRUE)
  dv_flat_2_s_w3 <- dv_flat_2_w3/sd(dv_flat_2_w1[Z=="control"],na.rm = TRUE)
  dv_flat_3_s_w3 <- dv_flat_3_w3/sd(dv_flat_3_w1[Z=="control"],na.rm = TRUE)
  dv_flat_4_s_w3 <- dv_flat_4_w3/sd(dv_flat_4_w1[Z=="control"],na.rm = TRUE)
  
  dv_mj_1_s_w3 <- dv_mj_1_w3/sd(dv_mj_1_w2[Z_distract==0],na.rm = TRUE)
  dv_mj_2_s_w3 <- dv_mj_2_w3/sd(dv_mj_2_w2[Z_distract==0],na.rm = TRUE)
  dv_mj_3_s_w3 <- dv_mj_3_w3/sd(dv_mj_3_w2[Z_distract==0],na.rm = TRUE)
  dv_mj_4_s_w3 <- dv_mj_4_w3/sd(dv_mj_4_w2[Z_distract==0],na.rm = TRUE)
  
  dv_climate_main_s_w3 <- dv_climate_1_s_w3
  dv_wall_main_s_w3 <- dv_wall_3_s_w3
  dv_amtrak_main_s_w3 <- dv_amtrak_1_s_w3
  dv_vets_main_s_w3 <- dv_vets_2_s_w3
  dv_flat_main_s_w3 <- dv_flat_1_s_w3
  
})

fa_climate <- psych::principal(filter(study_1_mturk, Z == "control") %>% select(dv_climate_1_s_w1, dv_climate_2_s_w1, dv_climate_3_s_w1, dv_climate_4_s_w1, dv_climate_5_s_w1), nfactors = 2, rotate = "varimax")
fa_wall <- psych::principal(filter(study_1_mturk, Z == "control") %>% select(dv_wall_1_s_w1, dv_wall_2_s_w1, dv_wall_3_s_w1, dv_wall_4_s_w1), nfactors = 2, rotate = "varimax")
fa_amtrak <- psych::principal(filter(study_1_mturk, Z == "control") %>% select(dv_amtrak_1_s_w1, dv_amtrak_2_s_w1, dv_amtrak_3_s_w1, dv_amtrak_4_s_w1), nfactors = 2, rotate = "varimax")
fa_vets <- psych::principal(filter(study_1_mturk, Z == "control") %>% select(dv_vets_1_s_w1, dv_vets_2_s_w1, dv_vets_3_s_w1, dv_vets_4_s_w1), nfactors = 2, rotate = "varimax")
fa_flat <- psych::principal(filter(study_1_mturk, Z == "control") %>% select(dv_flat_1_s_w1, dv_flat_2_s_w1, dv_flat_3_s_w1, dv_flat_4_s_w1), nfactors = 2, rotate = "varimax")
fa_mj <- psych::principal(filter(study_1_mturk, Z_distract == 0) %>% select(dv_mj_1_s_w2, dv_mj_2_s_w2, dv_mj_3_s_w2, dv_mj_4_s_w2), nfactors = 2, rotate = "varimax")

study_1_mturk <- within(study_1_mturk,{
  dv_climate_scale_w1 <- predict(fa_climate, data = study_1_mturk %>% select(dv_climate_1_s_w1, dv_climate_2_s_w1, dv_climate_3_s_w1, dv_climate_4_s_w1, dv_climate_5_s_w1))[,1]
  dv_wall_scale_w1 <- predict(fa_wall, data = study_1_mturk %>% select(dv_wall_1_s_w1, dv_wall_2_s_w1, dv_wall_3_s_w1, dv_wall_4_s_w1))[,1]
  dv_amtrak_scale_w1 <- predict(fa_amtrak, data = study_1_mturk %>% select(dv_amtrak_1_s_w1, dv_amtrak_2_s_w1, dv_amtrak_3_s_w1, dv_amtrak_4_s_w1))[,1]
  dv_vets_scale_w1 <- predict(fa_vets, data = study_1_mturk %>% select(dv_vets_1_s_w1, dv_vets_2_s_w1, dv_vets_3_s_w1, dv_vets_4_s_w1))[,1]
  dv_flat_scale_w1 <- predict(fa_flat, data = study_1_mturk %>% select(dv_flat_1_s_w1, dv_flat_2_s_w1, dv_flat_3_s_w1, dv_flat_4_s_w1))[,1]

  dv_climate_scale_w2 <- predict(fa_climate, data = study_1_mturk %>% select(dv_climate_1_s_w2, dv_climate_2_s_w2, dv_climate_3_s_w2, dv_climate_4_s_w2, dv_climate_5_s_w2))[,1]
  dv_wall_scale_w2 <- predict(fa_wall, data = study_1_mturk %>% select(dv_wall_1_s_w2, dv_wall_2_s_w2, dv_wall_3_s_w2, dv_wall_4_s_w2))[,1]
  dv_amtrak_scale_w2 <- predict(fa_amtrak, data = study_1_mturk %>% select(dv_amtrak_1_s_w2, dv_amtrak_2_s_w2, dv_amtrak_3_s_w2, dv_amtrak_4_s_w2))[,1]
  dv_vets_scale_w2 <- predict(fa_vets, data = study_1_mturk %>% select(dv_vets_1_s_w2, dv_vets_2_s_w2, dv_vets_3_s_w2, dv_vets_4_s_w2))[,1]
  dv_flat_scale_w2 <- predict(fa_flat, data = study_1_mturk %>% select(dv_flat_1_s_w2, dv_flat_2_s_w2, dv_flat_3_s_w2, dv_flat_4_s_w2))[,1]
  dv_mj_scale_w2 <- predict(fa_mj, data = study_1_mturk %>% select(dv_mj_1_s_w2, dv_mj_2_s_w2, dv_mj_3_s_w2, dv_mj_4_s_w2))[,1]
  
  dv_climate_scale_w3 <- predict(fa_climate, data = study_1_mturk %>% select(dv_climate_1_s_w3, dv_climate_2_s_w3, dv_climate_3_s_w3, dv_climate_4_s_w3, dv_climate_5_s_w3))[,1]
  dv_wall_scale_w3 <- predict(fa_wall, data = study_1_mturk %>% select(dv_wall_1_s_w3, dv_wall_2_s_w3, dv_wall_3_s_w3, dv_wall_4_s_w3))[,1]
  dv_amtrak_scale_w3 <- predict(fa_amtrak, data = study_1_mturk %>% select(dv_amtrak_1_s_w3, dv_amtrak_2_s_w3, dv_amtrak_3_s_w3, dv_amtrak_4_s_w3))[,1]
  dv_vets_scale_w3 <- predict(fa_vets, data = study_1_mturk %>% select(dv_vets_1_s_w3, dv_vets_2_s_w3, dv_vets_3_s_w3, dv_vets_4_s_w3))[,1]
  dv_flat_scale_w3 <- predict(fa_flat, data = study_1_mturk %>% select(dv_flat_1_s_w3, dv_flat_2_s_w3, dv_flat_3_s_w3, dv_flat_4_s_w3))[,1]
  dv_mj_scale_w3 <- predict(fa_mj, data = study_1_mturk %>% select(dv_mj_1_s_w3, dv_mj_2_s_w3, dv_mj_3_s_w3, dv_mj_4_s_w3))[,1]
  
  })


medians_df <- 
  study_1_mturk %>%
  summarise(climate_median = median(dv_climate_scale_w1[Z=="control"], na.rm = TRUE),
            wall_median = median(dv_wall_scale_w1[Z=="control"], na.rm = TRUE),
            amtrak_median = median(dv_amtrak_scale_w1[Z=="control"], na.rm = TRUE),
            vets_median = median(dv_vets_scale_w1[Z=="control"], na.rm = TRUE),
            flat_median = median(dv_flat_scale_w1[Z=="control"], na.rm = TRUE),
            mj_median = median(dv_mj_scale_w2[Z_distract == 0], na.rm = TRUE),
            
            climate_25 = quantile(dv_climate_scale_w1[Z=="control"], probs = 0.25, na.rm = TRUE),
            wall_25 = quantile(dv_wall_scale_w1[Z=="control"], probs = 0.25, na.rm = TRUE),
            amtrak_25 = quantile(dv_amtrak_scale_w1[Z=="control"], probs = 0.25, na.rm = TRUE),
            vets_25 = quantile(dv_vets_scale_w1[Z=="control"], probs = 0.25, na.rm = TRUE),
            flat_25 = quantile(dv_flat_scale_w1[Z=="control"], probs = 0.25, na.rm = TRUE),
            mj_25 = quantile(dv_mj_scale_w2[Z_distract == 0], probs = 0.25, na.rm = TRUE),
            
            climate_75 = quantile(dv_climate_scale_w1[Z=="control"], probs = 0.75, na.rm = TRUE),
            wall_75 = quantile(dv_wall_scale_w1[Z=="control"], probs = 0.75, na.rm = TRUE),
            amtrak_75 = quantile(dv_amtrak_scale_w1[Z=="control"], probs = 0.75, na.rm = TRUE),
            vets_75 = quantile(dv_vets_scale_w1[Z=="control"], probs = 0.75, na.rm = TRUE),
            flat_75 = quantile(dv_flat_scale_w1[Z=="control"], probs = 0.75, na.rm = TRUE),
            mj_75 = quantile(dv_mj_scale_w2[Z_distract == 0], probs = 0.75, na.rm = TRUE)
            )


study_1_mturk <- within(study_1_mturk,{
  dv_climate_agree_w1 <- as.numeric(dv_climate_scale_w1 >= medians_df$climate_median)
  dv_wall_agree_w1 <- as.numeric(dv_wall_scale_w1 >= medians_df$wall_median)
  dv_amtrak_agree_w1 <- as.numeric(dv_amtrak_scale_w1 >= medians_df$amtrak_median)
  dv_vets_agree_w1 <- as.numeric(dv_vets_scale_w1 >= medians_df$vets_median)
  dv_flat_agree_w1 <- as.numeric(dv_flat_scale_w1 >= medians_df$flat_median)
  
  dv_climate_agree_w2 <- as.numeric(dv_climate_scale_w2 >= medians_df$climate_median)
  dv_wall_agree_w2 <- as.numeric(dv_wall_scale_w2 >= medians_df$wall_median)
  dv_amtrak_agree_w2 <- as.numeric(dv_amtrak_scale_w2 >= medians_df$amtrak_median)
  dv_vets_agree_w2 <- as.numeric(dv_vets_scale_w2 >= medians_df$vets_median)
  dv_flat_agree_w2 <- as.numeric(dv_flat_scale_w2 >= medians_df$flat_median)
  dv_mj_agree_w2 <- as.numeric(dv_mj_scale_w2 >= medians_df$mj_median)
  
  dv_climate_agree_w3 <- as.numeric(dv_climate_scale_w3 >= medians_df$climate_median)
  dv_wall_agree_w3 <- as.numeric(dv_wall_scale_w3 >= medians_df$wall_median)
  dv_amtrak_agree_w3 <- as.numeric(dv_amtrak_scale_w3 >= medians_df$amtrak_median)
  dv_vets_agree_w3 <- as.numeric(dv_vets_scale_w3 >= medians_df$vets_median)
  dv_flat_agree_w3 <- as.numeric(dv_flat_scale_w3 >= medians_df$flat_median)
  dv_mj_agree_w3 <- as.numeric(dv_mj_scale_w3 >= medians_df$mj_median)
  
  dv_climate_agree_25_w1 <- as.numeric(dv_climate_scale_w1 >= medians_df$climate_25)
  dv_wall_agree_25_w1 <- as.numeric(dv_wall_scale_w1 >= medians_df$wall_25)
  dv_amtrak_agree_25_w1 <- as.numeric(dv_amtrak_scale_w1 >= medians_df$amtrak_25)
  dv_vets_agree_25_w1 <- as.numeric(dv_vets_scale_w1 >= medians_df$vets_25)
  dv_flat_agree_25_w1 <- as.numeric(dv_flat_scale_w1 >= medians_df$flat_25)
  dv_climate_agree_25_w2 <- as.numeric(dv_climate_scale_w2 >= medians_df$climate_25)
  dv_wall_agree_25_w2 <- as.numeric(dv_wall_scale_w2 >= medians_df$wall_25)
  dv_amtrak_agree_25_w2 <- as.numeric(dv_amtrak_scale_w2 >= medians_df$amtrak_25)
  dv_vets_agree_25_w2 <- as.numeric(dv_vets_scale_w2 >= medians_df$vets_25)
  dv_flat_agree_25_w2 <- as.numeric(dv_flat_scale_w2 >= medians_df$flat_25)
  dv_mj_agree_25_w2 <- as.numeric(dv_mj_scale_w2 >= medians_df$mj_25)
  dv_climate_agree_25_w3 <- as.numeric(dv_climate_scale_w3 >= medians_df$climate_25)
  dv_wall_agree_25_w3 <- as.numeric(dv_wall_scale_w3 >= medians_df$wall_25)
  dv_amtrak_agree_25_w3 <- as.numeric(dv_amtrak_scale_w3 >= medians_df$amtrak_25)
  dv_vets_agree_25_w3 <- as.numeric(dv_vets_scale_w3 >= medians_df$vets_25)
  dv_flat_agree_25_w3 <- as.numeric(dv_flat_scale_w3 >= medians_df$flat_25)
  dv_mj_agree_25_w3 <- as.numeric(dv_mj_scale_w3 >= medians_df$mj_25)
  
  dv_climate_agree_75_w1 <- as.numeric(dv_climate_scale_w1 >= medians_df$climate_75)
  dv_wall_agree_75_w1 <- as.numeric(dv_wall_scale_w1 >= medians_df$wall_75)
  dv_amtrak_agree_75_w1 <- as.numeric(dv_amtrak_scale_w1 >= medians_df$amtrak_75)
  dv_vets_agree_75_w1 <- as.numeric(dv_vets_scale_w1 >= medians_df$vets_75)
  dv_flat_agree_75_w1 <- as.numeric(dv_flat_scale_w1 >= medians_df$flat_75)
  dv_climate_agree_75_w2 <- as.numeric(dv_climate_scale_w2 >= medians_df$climate_75)
  dv_wall_agree_75_w2 <- as.numeric(dv_wall_scale_w2 >= medians_df$wall_75)
  dv_amtrak_agree_75_w2 <- as.numeric(dv_amtrak_scale_w2 >= medians_df$amtrak_75)
  dv_vets_agree_75_w2 <- as.numeric(dv_vets_scale_w2 >= medians_df$vets_75)
  dv_flat_agree_75_w2 <- as.numeric(dv_flat_scale_w2 >= medians_df$flat_75)
  dv_mj_agree_75_w2 <- as.numeric(dv_mj_scale_w2 >= medians_df$mj_75)
  dv_climate_agree_75_w3 <- as.numeric(dv_climate_scale_w3 >= medians_df$climate_75)
  dv_wall_agree_75_w3 <- as.numeric(dv_wall_scale_w3 >= medians_df$wall_75)
  dv_amtrak_agree_75_w3 <- as.numeric(dv_amtrak_scale_w3 >= medians_df$amtrak_75)
  dv_vets_agree_75_w3 <- as.numeric(dv_vets_scale_w3 >= medians_df$vets_75)
  dv_flat_agree_75_w3 <- as.numeric(dv_flat_scale_w3 >= medians_df$flat_75)
  dv_mj_agree_75_w3 <- as.numeric(dv_mj_scale_w3 >= medians_df$mj_75)
  
  
})



study_1_mturk <- within(study_1_mturk,{
  responded_w1 <- !is.na(dv_climate_scale_w1) & !is.na(dv_wall_scale_w1) & !is.na(dv_amtrak_scale_w1) & !is.na(dv_vets_scale_w1) & !is.na(dv_flat_scale_w1)
  responded_w2 <- !is.na(dv_climate_scale_w2) & !is.na(dv_wall_scale_w2) & !is.na(dv_amtrak_scale_w2) & !is.na(dv_vets_scale_w2) & !is.na(dv_flat_scale_w2)
  responded_w3 <- !is.na(dv_climate_scale_w3) & !is.na(dv_wall_scale_w3) & !is.na(dv_amtrak_scale_w3) & !is.na(dv_vets_scale_w3) & !is.na(dv_flat_scale_w3)
})

with(study_1_mturk, table(responded_w1))
with(study_1_mturk, table(responded_w2))
with(study_1_mturk, table(responded_w3))

study_1_mturk <- filter(study_1_mturk, responded_w1)

save(fa_climate,
     fa_wall,
     fa_amtrak,
     fa_vets,
     fa_flat,
     fa_mj,
     medians_df,
     file = "data/raw/factor_analysis_models.rdata")

write_rds(study_1_mturk, path =  "data/clean/study_1_mturk_cleaned.rds")


